function [B1 B2] = func_B_k(nz, nx, CONSTS)

    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    nxy = func_nxy_from_nz(nz, CONSTS);
    ny  = func_ny_from_nxy(nz, nx, CONSTS);
    [n1 n2] = func_nk_from_nz(nz, CONSTS);
    [J1 J2 J1_1 J2_1] = func_J_k(nxy, nx, CONSTS);

    B1 = -1i*(eta/eps_a)*J2_1.*(n2.*ny./(nxy.^2)) + (eta/eps_a)*J2.*(nx.*nz./(nxy.^2)) + ...
          n2.*((eps_a-nx.^2)./(eps_a*nxy.^2));
    B2 =  1i*(eta/eps_a)*J1_1.*(n1.*ny./(nxy.^2)) - (eta/eps_a)*J1.*(nx.*nz./(nxy.^2)) - ...
          n1.*((eps_a-nx.^2)./(eps_a*nxy.^2));
    
end
